Heterogeneous biological membranes regulate protein partitioning via fluctuating diffusivity

Abstract Cell membranes phase separate into ordered Lo and disordered Ld domains depending on their compositions. This membrane compartmentalization is heterogeneous and regulates the localization of specific proteins related to cell signaling and trafficking. However, it is unclear how the heterogeneity of the membranes affects the diffusion and localization of proteins in Lo and Ld domains. Here, using Langevin dynamics simulations coupled with the phase-field (LDPF) method, we investigate several tens of milliseconds-scale diffusion and localization of proteins in heterogeneous biological membrane models showing phase separation into Lo and Ld domains. The diffusivity of proteins exhibits temporal fluctuations depending on the field composition. Increases in molecular concentrations and domain preference of the molecule induce subdiffusive behavior due to molecular collisions by crowding and confinement effects, respectively. Moreover, we quantitatively demonstrate that the protein partitioning into the Lo domain is determined by the difference in molecular diffusivity between domains, molecular preference of domain, and molecular concentration. These results pave the way for understanding how biological reactions caused by molecular partitioning may be controlled in heterogeneous media. Moreover, the methodology proposed here is applicable not only to biological membrane systems but also to the study of diffusion and localization phenomena of molecules in various heterogeneous systems.

Biological membranes are composed of various kinds of proteins and lipids. Differences in the molecular composition relate to rich patterns of phase separation (1)(2)(3)(4)(5). Mixtures of saturated and unsaturated lipids generally cause phase separation into liquid-ordered (L o ) and liquid-disordered (L d ) domains (6)(7)(8). Specifically, the L d domain is rich in unsaturated lipids and of high fluidity, while the L o domain is rich in saturated lipids and of low fluidity. L o domains, enriched in sphingolipids and cholesterol, are often referred to as lipid rafts (9), and are thought to play a crucial role in a variety of cellular processes such as cell signaling and trafficking. Lipid rafts are generally considered to be small, heterogeneous, and highly dynamic domains of several tens of nanometers size with estimated life time 0.1-10 2 s (10)(11)(12)(13). The coexistence of L o and L d domains has been observed in synthetic model membranes under external stimuli or specific thermodynamic conditions. Direct imaging of cell-derived plasma membranes of giant plasma membrane vesicles has also confirmed the presence of nanodomains (14)(15)(16)(17)(18)(19)(20) at or near physiological temperature. Although there has been a longstanding debate regarding the nature and biological role(s) of these domains in living cells (21,22), a large number of recent studies have provided evidence for the coexistence of these domains in intact cells (23)(24)(25)(26)(27)(28)(29). These studies investigated the recruitment and exclusion of various probes associated with clustered proteins within cell membranes, and demonstrated that the concentration of probes in clusters reflects the partitioning observed in phaseseparated domains. The L o domains are formed not only by lipids but also by protein-lipid complexes, where the detailed properties, such as size, lifetime, and stability, depend on their composition and interaction with scaffolding proteins (21,22).
In terms of lateral diffusion of membrane proteins, phase separation may be considered as presenting an inhomogeneous field in which the protein molecules diffuse. The diffusivities of molecules in such inhomogeneous fields are known to be nonuniform in time and space (30,31). Experimental techniques, such as stimulated emission depletion (STED) microscopy combined with fluorescence correlation spectroscopy and single-particle tracking, have revealed dynamically heterogeneous motion of proteins in biological membranes (12,(32)(33)(34)(35)(36)(37)(38)(39)(40). Particularly, the local diffusivity of tracers fluctuates significantly with time due to the influence of heterogeneity in the field, e.g. intermittent trapping in domains (12,32), transient interactions with partners (34,38), or slow-active remodeling of the underlying cortical actin network (35,37). However, due to the difficulty of simultaneous measurement of molecular motion and field heterogeneity, the precise effects of membrane heterogeneity on molecular diffusion and distribution remain obscure. Although many theoretical models on molecular diffusion with fluctuating diffusivity have been proposed to explain the characteristics of non-Gaussian behavior and anomalous diffusion (41)(42)(43)(44)(45)(46)(47)(48)(49)(50)(51)(52)(53)(54)(55), it is important to understand the origin of fluctuations at the molecular level, specifically how phase separation, modeled as an inhomogeneous field, affects protein diffusivity and promotes protein crowding, or how molecular crowding induces phase separation and expands nanoscale domains. This understanding will help to clarify the role of proteinlipid and protein-protein interactions in the signaling process.
Molecular dynamics (MD) simulations have provided molecular details on protein diffusion in biological systems (56)(57)(58)(59)(60)(61)(62), and revealed temporal fluctuating of the protein diffusivity due to protein-protein and protein-lipid interaction (60,61). However, it remains a challenge for simulations to directly inform molecular dynamics on a spatiotemporal scale comparable to experiments. Here, using a mesoscale simulation technique, we unveil diffusion properties and distributions of molecules in heterogeneous biological membrane models. This coarse-grained level, combining Langevin dynamics simulations and phase-field methods, captures the motion of individual molecules in heterogeneous membranes at several tens of milliseconds timescales. We show the existence of fluctuating diffusivity and a distribution of molecules in heterogeneous membranes depending on various parameters such as heterogeneity of fields, molecular concentrations, and domain preference of molecules. This coarse-grained approach allows us to disentangle the effects of individual parameters on the observed protein motion, e.g. the diffusivity difference between the two membrane phases, the area density covered by proteins, or the protein affinity to a specific membrane domain. These results will be important to inform future experiments in real membrane systems in which some effects may be obscured by the complexity of the system.

Fluctuating diffusivity of an isolated molecule in heterogeneous biological membrane models
In our analysis, we focused on three distinct phase-separated heterogeneous biological membrane models described in previous studies (63). The phase separation process is measured in terms of the field c(r, t), the deviation of the local composition from the critical composition (see Materials and Methods section for the simulation details). The ordered (c < 0) and disordered (c > 0) phases denote the raft (L o ) and nonraft (L d ) domains, respectively. The distribution c(r, t) can be obtained by solving the reactiondiffusion equation. The specific model choices for the parameters induce clear phase separation and represent lipid raft formation; (model 2) interface pinning by immobile membrane proteins (64,65), (model 4) immiscible lipid systems, and (model 5) coupling to lipid reservoir (1, 66, 67) (see Fig. 1A). Considering the free energy term F, the phase separation is classified as "Mixed", "Nucleation", and "Spinodal Decomposition" (see Fig. 1B). Since the L o and L d domains have different compositions, diffusion coefficients of the biomolecules are different (12,13). To describe the diffusion of target protein molecules in such heterogeneous media, we considered the Langevin equation with fluctuating diffusivity, where r(t) is the position of a diffusing molecule at time t, and w(t) is white Gaussian noise with 〈w(t)〉 = 0. The diffusion coefficient D(r(t), t) varies depending on the field composition, is the normalized order parameter field (0 < c < 1) (see Fig. 1C and D for a sample trajectory). For a single molecular system, c b = 1 and D 0 = 1 were used in each model, i.e. D(r(t), t) fluctuates in the range of 1-2. The simulation time step dt = 0.001 and D 0 = 1 correspond to the physical quantities of 1 ns and 1 μm 2 /s, respectively. In the simulations, the system size L corresponds to 256 nm with periodic boundary conditions. Simulations were performed for 10 7 steps corresponding to 10 ms and analyzed after 10 6 steps (1 ms) of reaching equilibrium. Because the lifetime of the raft domain is 0.1-10 2 s (10-13), we here fixed the field variation and focused on timescales shorter than the field variation. This can allow us to evaluate the effect of spatial heterogeneity on the molecular dynamics, specifically elucidating how the (preexisting) raft domains affects the behavior of other molecules. First, we calculated the time-averaged mean squared displacement (TAMSD) as a quantity that characterizes the global diffusivity (see Fig. 2A), where Δ is a lag time and t is a measurement time. Individual TAMSDs increase linearly and show some amplitude scatter.
The probability density function (PDF) of TAMSDs at Δ = 10 −2 ms is found to have a distribution with roughly two peaks. This scatter is considered to be an effect of the inhomogeneity of the concentration distribution in the field. In order to quantitatively evaluate the effect of different patterns of heterogeneity on the diffusivity fluctuations, the relative standard deviation (RSD) of the TAMSDs was analyzed, It is known that RSD decays as t −0.5 in ergodic diffusion, e.g. Brownian motion. In the case of nonergodic diffusion processes, e.g. the continuous-time random walk (68,69), the RSD converges to a nonzero value for all Δ ≪ t as t → ∞. In fluctuating diffusivity models (43)(44)(45)70), the RSD exhibits a crossover from a plateau to a t −0.5 decay with a long crossover time. Here, the RSD shows a plateau in the time region t ∼ 10 −4 -10 −1 ms (see Fig. 2B), which implies that the instantaneous diffusivity fluctuates intrinsically on the corresponding timescale. Fluctuations of the diffusivity are negligible at the short and long timescales, where the RSD decays with t −0.5 . The short timescale depends on the initial diffusivity D(t = 0), while the long timescale relates to the relaxation time of the effective diffusivity. In a fluctuating diffusivity model where diffusivity dichotomously fluctuates between fast and slow states (44,45), the magnitude of the RSD depends on the difference in diffusion coefficients between the two states and the mean residence time of states. The magnitudes of the RSD of models 4 and 5 are higher than that of model 2.
To clarify the origin of the difference in RSDs, Fig. 2C shows the PDFs of the c for each model. The PDFs of models 2, 4, and 5 have two peaks and result in large diffusivity differences between L o and L d domains. Fig. 2D shows the PDFs of the residence times of the molecules in the L o and L d domains for each model. The residence times exhibit a power-law distribution with an exponential cutoff P(t) ∝ t −β exp ( − t/τ). The power-law exponents for each model are almost the same, β ≈ −1.5. The cutoff in the residence time relates to the relaxation time in the RSD at which the crossover from the plateau to the t −0.5 decay occurs. Longer residence times of the molecule in each domain translate into longer relaxation times of the RSD. Note that the first passage time (FPT) distribution of one-dimensional Brownian motion, starting from the origin at 0 and passing a certain point x, is given by the distribution , that is proportional to t −1.5 (t → ∞) (71), where D is the diffusion coefficient. When considering a finite-sized domain, the distribution t −1.5 has an exponential cutoff depending on the two-dimensional domain size. The general shape of the FPT distribution is similar for many scenarios (72).
We confirmed that slow variation of the concentration field affects little on the fluctuation of the diffusivity (see online supplementary Figs. S1 and S2). Since the timescale of the varying field is much longer than the simulation times, the domain boundaries change slightly in equilibrium states. In systems where the field varies faster than the timescale that particles move through the regions, a time-varying field may have a significant effect on the degree of the fluctuating diffusivity. In addition, we note that RSDs do not depend on the field patterns (see online supplementary Fig. S3).

Clustering effect of molecules on the fluctuating diffusivity in heterogeneous membranes
Cell membranes are crowded with a variety of proteins occupying 30-50% of the membrane area (73). In previous studies, a concentration dependency of protein subdiffusion, 〈δr 2 (Δ; t)〉 ∝ Δ α with α < 1, was observed in biological membranes (58,60). Switching off the protein-protein interactions changes the subdiffusive behavior (α = 0.84) to normal diffusion (α = 1.0) (34), dynamical correlations in the motions due to frequent molecular collisions may enhance subdiffusive motion (73). In any finite system, the subdiffusive regime will ultimately cross over beyond some correlation time, see, e.g. Ref. (74). To explore the effect of membrane crowding, we evaluate the diffusivity of molecules in molecular crowded systems with N p = 64, 128, 256, 512, 1024, and 2048 particles corresponding to an area occupancy of 1.4, 3.5, 7.8, 16.6, 34.0, and 59.1% of the L o domains, respectively. In the following, we mainly focus on model 5 (results for other models are shown in online supplementary Fig. S4). In this membrane state, the separation of L o and L d domains is most pronounced and thus best accessible in experiments. Fig. 3A shows the aggregation of molecules with different area occupancy (see online supplementary Movie S1). Even in the absence of molecular field preference, we find that as N p increases, molecules tend to aggregate in the L o domain, where the diffusion coefficient of molecules is smaller than in the L d domain. This aggregation affects the diffusive behavior of molecules. Ensemble-averaged TAMSDs become smaller and exhibit subdiffusion when the area occupancy increases (see Fig. 3B). The power-law exponent of the TAMSD decreases from α = 1.0 to 0.85, depending on the molecular concentration, up to a timescale of ∼0.1 ms. This trend is similar to that of MD simulations reporting transient subdiffusion of proteins in a molecular concentrationdependent manner (58,60). Coarse-grained MD simulation for 0.1 ms (60) showed that subdiffusive motion of proteins in a noncrowded membrane changes to Brownian motion at Δ > 10 ns attributed to the viscoelasticity of lipids (75)(76)(77)(78), while in a crowded membrane significant subdiffusive regimes α ∼ 0.8-0.9 extends until tens of microseconds (>0.01 ms is not conclusive due to the limited simulation time). Similar transient subdiffusion induced by molecular crowding was also observed for hard-core particles (79).
Here, we discuss the effect of the concentration and the interaction strength of the tracer molecules. Molecular concentration has little effect on the magnitude of RSD, and the relaxation time increases slightly with increasing molecular concentration. The molecular concentration differences have little effect on the fluctuation of the diffusivity (see Fig. 3B). Moreover, we evaluated the effect of the interaction strength ϵ between molecules (see Eq. 9 in Materials and Methods section). When the molecular interaction becomes stronger, the magnitude of the TAMSD becomes small and α decreases (see Fig. 3C). An increase in the interaction strength has little effect on the RSD.

Preference of the domain affects the diffusivity in heterogeneous membranes
It is known that the differences in lipid compositions in L o and L d domains generate a preferential partitioning of membrane proteins in either domain. The protein domain preference, especially of transmembrane proteins, is determined by its palmitoylation, hydrophobic length, and surface area of its transmembrane region (80,81).
In our simulations, the preference was modeled using a reflective wall at the boundary between L o and L d domains (see details in Materials and Methods section). We evaluate the effect of preference of the L o domain (L o χ) (Fig. 4A) or the L d domain (L d χ) (Fig. 4B) on the diffusive dynamics, where χ is the degree of the domain preference. As shown in Fig. 4A, molecules are localized more in the L o domain with strong L o domain preference. According to an increase of χ, the TAMSD decreases, and the molecules exhibit more pronounced subdiffusion with smaller anomalous exponents α = 0.8-1.0. In the case of L d domain preference, molecules are localized more in the L d domain, and the TAMSD increases with higher χ (Fig. 4B). Molecules exhibit subdiffusion with anomalous exponents α = 0.9-1.0. Note that the crossover of α < 1 to normal diffusion is not observed for larger L o χ in the studied time window (Fig. 4A). This means that the caging effect of molecules in narrow regions strongly influences anomalous diffusion, significantly more than the crowding effect with high concentrations (timescale of ∼0.1 ms in Fig. 3).
The magnitude of the RSDs for both L o χ and L d χ becomes smaller upon increase of χ (Fig. 4). This is thought to be due to the fact that high χ increases the confinement of molecules to a preferable domain, which leads to a decrease in the fluctuation of diffusivity. Moreover, the area of the L d domain is larger than that of the L o domain in model 5. The residence time of the molecule increases with growing domain area, and the diffusivity of the molecule remains the same, resulting in a smaller RSD value. Note that in model 2 and model 4 membranes, where the areas of L d and L o domains are the same (see Fig. 2C), the change in RSD when domain preference is changed is almost the same for L o χ and L d χ (see online supplementary Fig. S4).

Confinement of molecules to one domain due to membrane heterogeneity
A nanoscale domain in membranes increases local molecular concentrations and molecular collisions, which are relevant to biological reactions. To see this, the distribution of molecules in the heterogeneous membrane was analyzed. Fig. 5 shows ratios of molecules confined in the L o domain examined for each parameter, such as N p , ϵ, and domain preference. Randomly distributed particles at the initial time (t = 0) diffuse and start to enrich in the L o domain times of t = 0.01 to 0.1 ms. The confinement ratio changes like a sigmoidal curve and reaches a plateau (equilibrium) after 0.1 ms (see Fig. 5A and B). Although there is no preferential domain for molecules (L o 0 and L d 0), molecules are more likely to stay in the L o domain, where the diffusion coefficient is lower than in the other domain, and aggregate with surrounding molecules there. An increase in molecular concentration enhances the speed of the ratio increase and the equilibrated ratio because of high encounter rates at high concentrations (see Fig. 5A). The confinement into the L o domain is also enhanced by the interaction strength ϵ between molecules (see Fig. 5B). The strength of ϵ does not affect the speed of the ratio increase but increases the ratio at the plateau (t > 1 ms) as a high interaction strength stabilizes the cluster of molecules.
We now examine the effect of domain preference, L o χ or L d χ. According to an increase in the degree of preference χ, once the molecules enter the preferable domain, molecules cannot easily exit from the domain. Fig. 5C shows that an increase of χ of L o domain increases both the equilibrated ratio and the speed of the ratio. While an increase of χ of L d decrease the ratio of molecules in the L o domain with a crossover around χ ∼ 60 (see Fig. 5D). The preferential distribution of molecules to the domain of low diffusivity is inverted by the affinity strength between molecules and the domain of high diffusivity.

Discussion
Visualizing small and highly dynamic domains in cell membranes can be challenging due to the reduction in binary contrast of the heterogeneity when averaging over time. Despite prolonged debate (21,22), recent studies have provided compelling evidence for the coexistence of L o and L d domains in living cells (23)(24)(25)(26)(27)(28)(29). The mobility and aggregation of membrane proteins in cell membranes are closely linked to the local lipid order, i.e. phase separation is an organizing principle for membrane protein partitioning (27). Recent experiments using a broad range of fluorescent probes with various membrane anchors, which have different lipid domain preferences, show the existence of segregated domains selectively partitioning membrane proteins according to their affinity for the L o or L d domain (29). Some membrane proteins are enriched in the nanoscale region surrounding the clustered receptors. Although these studies on multicomponent systems provide us with insights into macroscopic biological regulation through heterogeneity in membranes, the underlying intricate mechanisms regulating the molecular dynamics in such complex systems have not been fully dissected. A theoretical understanding is crucial to gain insights into the precise factors that control molecular diffusion and localization within inhomogeneous fields. Such a quantitative model is also indispensable for data analysis in membrane systems with advanced assimilation methods based on prior training (82,83). In this study, we have used a well-defined in silico setup simulating molecules with fluctuating diffusivity in phase-separated fields with L o (low diffusivity) and L d (high diffusivity) domains. This coarse-grained approach allows us to disentangle the various effects conspiring in the complex observed motion. We showed that the degree of fluctuating diffusivity depends on the magnitude of the difference in molecular diffusivity between domains and the residence time in domains. Our results suggest that molecular localization within L o (low diffusivity) domains spontaneously occurs in heterogeneous membranes even when there is no domain preference, and subdiffusive behavior is observed due to molecular collision via molecular crowding in L o domains. Domain preference extends the timescale of the subdiffusive regime via molecular confinement into the preferential domains. The effect of heterogeneity on protein partitioning was also quantitatively evaluated. We demonstrated that the localization of molecules is determined by the difference in molecular diffusivity between domains, molecular preference of domain, and molecular concentration.
The aforementioned results were obtained under the condition of a fixed field variation to dissect the effect of preexisting raft domains on the molecular behavior. Membrane proteins possess the ability to modify their surrounding lipid environment, leading to the formation of functional protein-lipid complexes. Our approach could also be applicable in scenarios where protein-lipid interactions promote or alter the functional domains within the membrane. To model such effects, we conducted simulations under two distinct conditions: one representing a scenario where macroscopic phase separation does not occur spontaneously without proteins (near the miscibility critical point, Λ = −0.01) (Fig. 6A), and another representing a state where macroscopic phase separation occurs spontaneously without proteins (under the miscibility critical point Λ = −1) (Fig. 6C). Fig. 6A shows snapshots of the phase-separated field driven by diffusing proteins. Even at a condition where the macroscopic phase separation does not occur spontaneously, from a uniformly distributed state of the field and proteins, the formation of small-scale clusters of proteins leads to local L o domains. At this simulation condition, the normalized c has clear one distinct peak around ∼0.6 (L d ) and rudder point at c < 0.5 (L o ) (Fig. 6B). Note that the contrast between L d and L o and the domain size depend on the protein-protein and protein-lipid interaction strength.
At a condition under the critical miscibility temperature, macroscopic phase separation into L d and L o occurs with marked contrast (Fig. 6C). The PDF of c has two peaks around ∼0.8 (L d ) and ∼0.2 (L o ) (Fig. 6E). This relates to a scenario of the recruitment of additional proteins to the existing functional domains and their subsequent alteration of the domain configuration and function. Interestingly, the presence of diffusing molecules causes the distorted domain configuration, while the absence of molecules results in the formation of spherical domains (Fig. 6D). In addition, diffusing molecules cause the extension and fusion of the domains (Fig. 6F). Such recruitment can alter the thermodynamic stability of the membrane domains without a change in the lipid composition. MSD (Fig. 6G) and RSD (Fig. 6H) are consistent with the previous results (fixed field variation) that magnitude of the RSD depends on the difference in diffusion coefficients between the L o and L d domains.
In realistic biological membranes, lipid compositions vary significantly for different cell types. External ions and biomolecules further moderate the heterogeneity in the signaling and trafficking processes. These factors regulate the heterogeneity of the phase-separated membrane and formation of functional multiprotein units in membranes (84). Interaction with the underlying actin cytoskeleton also regulates condensation of the phase along the actin filament by pinning elements to a preferred phase (85). In suitable conditions, the field-dependent diffusive behavior of molecules is expected to regulate the search time of partners and reaction rates (86). Protein condensation on the phase-separated membrane surfaces is a key role in downstream signaling (87,88). A quantitative and qualitative elucidation of the nature of molecular behaviors in heterogeneous media is critical to understanding cellular behavior.
Our approach presented here is quite general and can be applied to fundamental questions on molecular dynamics in a variety of heterogeneous media in biology, soft matter, solidstate physics, etc. Our model could also be extended to more realistic biological membrane models including, e.g. dynamic modulation of protein domain preferences via phosphorylation by interaction with regulatory proteins (89), protein remodeling through conformational changes, complex inhomogeneous interaction between molecules and clusters, the partitioning by an actin filament mesh (37,90), alternation of membrane composition in signaling events, and the partitioning regulation of membrane signaling (91). Moreover, the parameters of this mesoscale simulation can be determined bottom-up from MD simulations, allowing comparison of mesoscopic molecular behavior at the intersection of simulations and experimental spatiotemporal scales. Concurrently, advanced single particle tracking studies provide massive new data on protein dynamics in membranes that can be scrutinized by our approach and advanced methods for data analysis (82,83). This could open a new direction to delineate the role of heterogeneity in the membranes with more complex multicomponent systems in a more physiological setting (92,93).

Simulation models
We used five models of phase separation in cell membranes as described in Ref. (63). The specific choices for the parameters represent lipid raft formation (model 1) by thermal fluctuations near the critical temperature, (model 2) by pinning of the interfacial composition of immobile membrane proteins (64,65), (model 3) in miscible or (model 4) immiscible lipid systems, and (model 5) by exchange with lipid reservoirs (1,66,67). These five models are expressed using a Cahn-Hilliard equation (63,94,95) for the order parameter field c(r, t), The first term on the right-hand side is the term for the lipid reservoir in model 5, where τ r is a parameter representing the relaxation time due to coupling to the lipid reservoir, and c r is the average compositions imposed by the lipid reservoir. The second term on the right-hand side is the modified Ginzburg-Landau free energy term in the usual Cahn-Hilliard equation, where M is the mobility and F is the free energy, where the parameter Λ > 0 (Λ < 0) is a relative temperature to the mean-field critical temperature T > T c (T < T c ). W is a parameter to control the line tension between the raft and nonraft phases. α is a parameter that explain the local reduction in the line tension due to immobile membrane proteins. The local concentration ρ(r) of N immobile membrane proteins in model2 can be expressed as η(r, t) in Eq. 4 denotes a Gaussian noise term (94), where H can be related to either the temperature of the system or the rate at which lipids are removed and added to the leaflet due to vesicular and nonvesicular lipid trafficking events, l denotes the recycling length over which spatial redistribution of lipids takes place (96), and ξ (q, t) is the Fourier transform of the white Gaussian noise with mean 0 and variance 1. Here, we used , W = 1, and the values of each parameter in each membrane model are shown in Table 1 (63). The ordered (c < 0) and disordered (c > 0) phases denote the raft (L o ) and nonraft (L d ) domains.
The system was simulated using the phase-field method under periodic boundary conditions with a grid point size of 256 × 256 (256 nm × 256 nm in physical dimensions). The lattice point width was set to Δx = Δy = 1 (1 nm for physical quantities). The time step was set to dt = 0.005 for dimensionless numbers, which corresponds to 10 −5 s for physical quantities. The number of simulation steps for each model is shown in Table 1. The initial c field was set to be homogeneous with a Gaussian distribution with mean 0 and variance 1.

Single-particle system
The diffusive particles in each membrane model are modeled by the Langevin equation 1 with fluctuating diffusivity. The diffusivity of the particle, D(r(t), t) = (c b + c(r(t), t))D 0 , fluctuates depending on the normalized order parameter field c(r(t), t) (0 < c < 1). We used c(r(t), t) in equilibrium after running simulations for each number of steps in Table 1. The single-particle simulations were performed 100 times with different initial coordinates of the particles for the same phase-separated field. The parameters c b = 1 and D 0 = 1 (1 μm 2 /s) were used in each model. Simulations were performed for 10 7 steps (10 ms) with dt = 10 −3 (1 ns), and the trajectories of the particles were analyzed after 10 6 steps (1 ms) of reaching equilibrium.

Multiparticle system
For multiparticle interactions, we performed simulations including particle-particle interactions, where k B T = 1, and Lennard-Jones potential was used, where l was the distance between two interacting particles, size of the particle σ was 3.0. The depth of the potential well ϵ was set as 0.5, 1.0, 2.0. The number of particles in the system N p was set to N p = 64, 128, 256, 512, 1, 024, 2, 048 to compare the effect of particle concentration on the diffusivity. For multiple particle systems, we used c b = 0.1 and D 0 = 1 in D(r(t), t) = (c b + c(r(t), t))D 0 . We used cdview (https://polymer.apphy.u-fukui.ac.jp/~koishi/ cdview.php) for visualization of the simulations.

Nanoscale phase separation model modified by diffusive particles
For a scenario where protein-lipid interaction induces nanoscale phase separation and creates functional domains in the membrane, we conducted a simulation in which diffusing particles change the field. For phase separation, we considered the following expression: where M = 1, g(r, t) is a term of short-ranged protein-lipid interaction, g(r, t) = α g (r ≤ σ), g(r, t) was considered at the positions of the diffusing particles.
Here, we set the intensity of field modification by proteins as α g = 0.5 and its relaxation length as r g = 2.

Domain preference of molecules
To implement the domain preference of the molecule, the energy barrier between the domains was reproduced by probabilistic reflection when a particle moves from one domain to another. We compared three patterns. One is that the particles can move freely between the L o and L d domains without being reflected (χ = 0). The other two are cases where the particles exhibit L o or L d preferences. L o χ means that the molecule is reflected at L o when moving from the L o domain to the L d domain at χ% probability of reflection and not reflected when moving in the opposite direction. L d χ is vice versa.

Supplementary material
Supplementary material is available at PNAS Nexus online.

Author contributions
E.Y. designed the research; K.S. performed the simulations; K.S. and E.Y. analyzed the data; the research emerged from discussions between all authors; K.S. and E.Y. mainly wrote the manuscript, and all authors partially contribute to write the manuscript.

Data availability
All relevant data are provided in the article and/or online supplementary material.